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geometry from spacecraft data: MHD simulation with guide field 
and antiparallel kinetic simulation 
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[i] When analyzing data from an array of spacecraft (such as Cluster or MMS) crossing a 
site of magnetic reconnection, it is desirable to be able to accurately determine the 
orientation of the reconnection site. If the reconnection is quasi-two dimensional, there are 
three key directions, the direction of maximum inhomogeneity (the direction across the 
reconnection site), the direction of the reconnecting component of the magnetic field, and 
the direction of rough invariance (the “out of plane” direction). Using simulated spacecraft 
observations of magnetic reconnection in the geomagnetic tail, we extend our previous 
tests of the direction-finding method developed by Shi et al. (2005) and the method to 
determine the structure velocity relative to the spacecraft Ktr- These methods require data 
from four proximate spacecraft. We add artificial noise and calibration errors to the 
simulation fields, and then use the perturbed gradient of the magnetic field B and perturbed 
time derivative dB/dt, as described by Denton et al. (20 1 0). Three new simulations are 
examined: a weakly three-dimensional, i.e., quasi-two-dimensional, MHD simulation 
without a guide field, a quasi-two-dimensional MHD simulation with a guide field, and a 
two-dimensional full dynamics kinetic simulation with inherent noise so that the apparent 
minimum gradient was not exactly zero, even without added artificial errors. We also 
examined variations of the spacecraft trajectory for the kinetic simulation. The accuracy of 
the directions found varied depending on the simulation and spacecraft trajectory, but all 
the directions could be found within about 10° for all cases. Various aspects of the method 
were examined, including how to choose averaging intervals and the best intervals for 
determining the directions and velocity. For the kinetic simulation, we also investigated in 
detail how the errors in the inferred gradient directions from the unmodified Shi et al. method 
(using the unperturbed gradient) depended on the amplitude of the calibration errors. For 
an accuracy of 3° for the maximum gradient direction, the calibration errors could be as 
large as 3% of reconnection magnetic field, while for the same accuracy for the minimum 
gradient direction, the calibration errors could only be as large as 0.03% of the reconnection 
magnetic field. These results suggest that the maximum gradient direction can normally be 
determined by the unmodified Shi et al. method, while the modified method or some other 
method must be used to accurately determine the minimum gradient direction. The structure 
velocity was found with magnitude accurate to 2% and direction accurate to within 5%. 
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1. Introduction 

[ 2 ] One of the major goals of magnetospheric spacecraft 
observations is to investigate spatial structures such as sites of 
magnetic reconnection. The starting point for such investi- 
gations is to determine the orientation of the structure relative 
to the spacecraft. While there has been considerable research 
on methods to determine the direction across a one dimen- 
sional discontinuity from spacecraft observations [ Sonnerup 
et al., 2006, and references therein; see also Mozer and 
Retind, 2007], there has been less research on methods to 
determine the orientation of two-dimensional structures. In 
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Figure 1 . Sketch of reconnection geometry in the magnetotail. The X point is at the origin of the X-Z (GSM) 
coordinate system. At large |Z|, the reconnecting magnetic field is oriented in the X direction. Magnetic flux is 
transported toward the X point from above and below with speed v in , and away from the X point to the right 
and left with speed v out . The width 8 of the diffusion region (rectangular box), is less than its length L, 
corresponding to greater gradient in the Z direction. The path of the centroid of the array of virtual spacecraft 
is schematically represented by the path S. (The spacecraft separation is small on this scale.) 


this paper, we will consider “quasi-two dimensional” recon- 
nection, as illustrated in Figure 1. If a structure is quasi-two 
dimensional, there will be a direction in which the variation 
is significantly less than that of the other directions. We 
call this the “invariant” or “out of plane” direction, recog- 
nizing however that in a real three-dimensional system, 
this invariance is only approximate. In the case of quasi- 
two dimensional reconnection the equilibrium (large-scale) 
gradients within the plane shown in the figure will be 
much larger than the gradients out of the plane. The direction 
of largest gradient in Figure 1 is across the current sheet, 
i.e., the Z direction (e z ) in the figure. The direction of inter- 
mediate gradient is along the reconnecting magnetic field, 
i.e., the X direction (e x ) in the figure. The direction of 
minimum gradient is the third (out-of-plane) direction, i.e., 
the Y direction (e y ). 

[ 3 ] Shi et al. [2005] developed a method to determine 
three directions corresponding to maximum, intermediate, 
and minimum values of the squared magnitude of the vector 
magnetic field gradient. They called this technique “Mini- 
mum Direction Derivative (or Difference)”, or MDD, anal- 
ysis. This method requires field observations from four 
closely spaced spacecraft in order to determine the gradient. 
A potential advantage of this method is that it yields the 
directions at each point in time, thus in principle allowing 
changes in orientation of a structure to be monitored. 
According to the method, the gradient of the magnetic field 
is first calculated, and expressed as the matrix M V8 , where 
M^ 8 = djB k , and d t is the spatial partial derivative in the f th 
direction. Then the symmetric matrix M G = M V8 • (M V8 ) T 
is formed, where “T” indicates the transpose of the matrix. 
The three eigenvalues of M°, A G _ max , A G _ int , and A G _ min , 
represent the maximum, intermediate, and minimum values 
of the squared directional derivative (gradient), with the 
eigenvectors <? G _ max , e G _ int , and <? G _ min indicating the 
corresponding directions. As Shi et al. [2005] explain, in 
order for the structure to be roughly two-dimensional, 
A G — min should be significantly less than the other two 
eigenvalues. If A G _ max » A G _ int , A G _ min , the structure is 
quasi-one-dimensional. If all the eigenvalues are 


comparable, the structure is fully three-dimensional, and the 
eigenvector directions may not be useful for our purposes. 
Shell et al. [2003, 2007a, 2007b] used a similar approach 
based on the gradient of the magnetic field direction h = BIB 
(see discussion by Denton et al. [2010]). Note that the Shi 
et al. [2005, 2006] method has been used recently to inves- 
tigate small-scale magnetic structures [Shi et al., 2009a] and 
the cusp boundary [Shi et al., 2009b]. 

[ 4 ] Assuming time stationarity (d/dt = 0) in the frame of 
the structure, Shi et al. [2006] went on to use dB/dt and Vi?, 
observed by the four spacecraft to determine the velocity of 
a structure relative to the spacecraft, V str = — V sc , where V sc 
is the velocity of the spacecraft array relative to the structure, 


dB 

dt 


dB 

"X f Vsc ' Vi? = 

at 


v str Vfi, 


(1) 


They called this technique “Spatio-temporal Difference” 
analysis, which they abbreviated as STD. 

[ 5 ] Denton et al. [2010] considered the effect of random 
noise and calibration errors on the Shi et al. [2005, 2006] 
methods. While the effect of noise errors could be elimi- 
nated by averaging in time, calibration errors were a more 
difficult problem to deal with because they systematically 
contaminate calculation of the gradient V /? [Denton et al., 
2010]. This is because the magnetic field measured by the 
spacecraft is the sum of the real magnetic field and the con- 
stant (with respect to time) calibration errors, and the gradient 
of the calibration errors is constant and nonzero, leading to 
systematic error in the apparent gradient. 

[6] In order to eliminate the effect of calibration errors, 
Denton et al. [2010] modified the Shi et al. [2005] method 
by using the perturbed matrix 

SM V8 = M v *-(m V8 ) 0 , (2) 

where (M v ' s } 0 is an average value of M V8 (=Vi?) calcu- 
lated within an interval typically centered on the time at 
which M V8 has its largest value. This is also the time of 
largest eigenvalue of M , since that eigenvalue is the square 
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of the amplitude of the largest value of V B [ Denton et al., 
2010]. They then used <5M v/i for the analysis in the 
same manner as Shi et al. [2005]. Subtracting off the 
average value of M v/i eliminates the effect of the calibration 
errors, since these errors lead to a constant gradient. Denton 
et al. then found the time average (<5M G ) of the time depen- 
dent matrix <5M G = <5M VS • (<5M V8 ) T , where <5M V8 = M V8 - 
(M vs ) 0 . The averaging interval used to calculate (bM G ) may 
be different from that used for (M V8 ) 0 . Experience to date 
indicates that it may be best to average (M V8 ) 0 near the 
center of the current sheet where the gradient is large, 
whereas the minimum gradient direction is sometimes best 
found away from the current sheet [Denton et al., 2010], 
Having averaged (<5M G ), we can then use the eigenvectors of 
(bM G ) to get the gradient directions. By use of shifting 
averages, slow changes in orientation can in principle still be 
monitored, although this aspect is not pursued here. 

[ 7 ] Denton et al. [2010] also modified (1) in order to 
eliminate the effect of calibration errors. Using the average 
of VB and the average of dB/dt (in the spacecraft frame) 
over the same time interval, V str can be found from 

=~Y sl rS(VB), (3) 

where b(dB/dt ) = dB/dt — ( dB/dt) q , and 6(VB) = VB — (VB) 0 . 

[8] Using simulated spacecraft observations from four 
virtual spacecraft flown through an MHD simulation of 
magnetic reconnection in the geomagnetic tail, Denton et al. 
[2010] tested the direction-finding method [Shi et al., 2005] 
and the method to determine the structure velocity [Shi et al., 
2006]. Employing the gradient calculated from the measured 
vector magnetic field without calibration errors and noise, 
they found that these quantities could be well determined. 
They also showed that these quantities could be well deter- 
mined even when calibration errors were added, provided 
that the perturbed quantities were used in the calculations as 
discussed above. They found that the angles were deter- 
mined within about 3° and the magnitude of the structure 
velocity within about 3%. 

[ 9 ] While the results of Denton et al. [2010] are encour- 
aging, there are limitations of the study that we address in this 
paper. First, the results were derived entirely from a single 
MHD simulation of reconnection in the magnetotail. We 
now want to see if the method works for other simulations. 
The simulation used by Denton et al. was for anti-parallel 
reconnection (no component of the equilibrium magnetic 
field in the out of plane e Y direction). In the present paper, we 
will test the method for another quasi-two dimensional, anti- 
parallel MHD reconnection simulation, but we will also test 
it using an MHD simulation with a guide field and a full 
dynamics kinetic simulation. 

2. Errors 

[ 10 ] As discussed by Denton et al. [2010], in real space- 
craft data there are two kinds of errors in magnetic field 
measurements, digitization (noise) errors that randomly vary 
with respect to time, and systematic calibration errors that 
are constant or very slowly evolving with respect to time. 
For each simulation, following Denton et al., we will add to 
the virtual spacecraft measurements 0.01 nT time varying 


noise errors and constant 0.1 nT calibration errors, both 
randomly chosen for each component of the magnetic field. 
(For calibration errors of this magnitude, the apparent value 
ofV-B (from the trace of the M^ 8 matrix) is not signifi- 
cantly affected, though it could be for larger calibration errors.) 
We consider this a worst-case scenario. For Cluster, it is 
sometimes possible to reduce the calibration errors in the 
spacecraft spin plane using measurements by the Electron 
Drift Instrument (EDI) [Georgescu et al., 2006]. For MMS, 
it should be possible to reduce the relative (between space- 
craft) calibration errors (which is crucial for calculating 
the gradient), taking advantage of time periods for which the 
magnetic field measured by each spacecraft should be the 
same because the field is nearly uniform and the spacecraft 
separation is small (R. Torbert, private communication, 
2010). Reduction of the errors will make the calculations 
more accurate. For the present paper, we will ignore errors 
associated with calculation of the gradients; basically, the 
spacecraft configuration needs to be roughly tetrahedral and 
the scale length for spatial variation needs to be significantly 
larger than the spacecraft separation in order for the gra- 
dients to be calculated accurately. 

3. Investigations 

[ 11 ] Denton et al. [2010] considered a snapshot of the 
MHD fields from a quasi-two-dimensional simulation of 
reconnection [Birn and Hesse, 2009]. On the large scale 
[Birn and Hesse, 2009, Figure 3], the structure was strongly 
three dimensional, but at the central value of the simulation 
(y = 0 [Birn and Hesse, 2009, Figure 3], where y is the cross 
tail (GSM Y) coordinate), the simulation was approximately 
two dimensional. That is, their B v was small compared to the 
lobe magnetic field and the cross-tail variation was weak 
compared to the variation in their z (downstream tail coor- 
dinate, GSM —X). Virtual (simulated) spacecraft were flown 
across this snapshot (with a path schematically represented 
by S in Figure 1) so that the only time dependence of the 
recorded data was due to the motion of the spacecraft. (The 
spacecraft velocity included components in all three direc- 
tions including the quasi-two-dimensional “out-of-plane” 
direction (y direction of Birn and Hesse).) This same method 
will be used in this paper; that is, real time dependence will 
not be considered. Therefore, the potential of the Shi et al. 
method to recover slow time evolution of directions and 
velocity will not be examined. Anomalous resistivity was 
used in the region around the X point which was traversed 
by the spacecraft. The reconnecting magnetic field (the X 
component) was 20 nT and the proton density was of order 
0.08 cnU 3 in the current sheet, and 0.05 cm~ 3 in the 
upstream (inflow) region. The plasma beta was small ~10~ 3 
to 10~ 2 in the inflow region, but large »1 at the current 
sheet. This simulation had approximately zero guide field 
(zero Y component of B in the midnight meridian plane). 
Here we will first consider data from two other MHD 
simulations with similar parameters. 

3.1. MHD Simulation With Zero Guide Field 

[ 12 ] Run 1 is very much like the simulation considered by 
Denton et al. [2010]; it also does not have a guide field. The 
main difference between the two simulations is that the 
length scale in the new simulation is smaller (200 km versus 
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t(s) 


Figure 2. The x (black), y (blue), and z (green) components 
(rotated in an arbitrary way relative to the natural simulation 
coordinates) of the magnetic field for all four virtual space- 
craft flown through the snapshot of the magnetic field from 
Run I. 


1000 km), so that the crossing of the current sheet is faster 
and there are consequently fewer data points during the 
crossing. Before analyzing the data, we add 0.01 nT random 
noise error (varying with time) and 0.1 nT systematic cali- 
bration errors (constant with respect to time) to each com- 
ponent of the magnetic field. The data were sampled at a 
time resolution of 1.0 s. Figure 2 shows the magnetic field 
measurements in (arbitrary) simulation coordinates for all 
four virtual spacecraft. The fact that the four curves for each 
component do not vary greatly (the blue curves, for instance, 
are almost exactly superposed) indicates that the separation 
of the spacecraft (100 km here) is small relative to the global 
spatial scales. To implement the Shi et al. [2005] method 
using the perturbed magnetic field developed by Denton 
et al. [2010], we first need to identify the region of the 
largest gradient, which occurs when the central current sheet 
is crossed at approximately t = 0. Figure 3 shows the three 
eigenvalues, A G - ma x (black), A G _ int (blue), and A G _ min 
(green) using the matrix M VB based on the gradient of the 
unperturbed magnetic field after doing a running average of 
the magnetic field values at each time using the values 
within ±5 s (1 1 measurements) in order to reduce the effects 
of noise. The peak in A G _ max occurs at about t = — 3 s. Using 
a procedure to be discussed shortly, we chose a time seg- 
ment between —41.5 s and 34.5 s for the purpose of calcu- 
lating (M v,s >o in order to get <5M VB using (2). 

[ 13 ] Figure 4a shows the eigenvalues using the perturbed 
gradient matrix <5M V ' S , again after doing a running average 
of the magnetic field values at each time using the values 
within ±5 s. As mentioned in the Introduction, averaging the 
data helps to eliminate problems related to random noise, 
whereas using the perturbed gradient <5M VB eliminates the 


problem of calibration errors. Note that the maximum gra- 
dient eigenvalue A G _ max (black) in Figure 4a is now much 
lower in the central region, because we have subtracted off 
the average gradient within this region. Nevertheless, this 
eigenvalue is still significantly higher than the other eigen- 
values except at two times where A G _ max plummets to small 
values; this is where the value of A G _ max based on the 
unperturbed gradient approaches the average value of the 
gradient within the central region. Figures 4b^ld show 
the direction cosines for the maximum gradient direction 
^G— nmx! the intermediate gradient direction e G _ int , and the 
minimum gradient direction e G _ m ; n . The gradient directions 
are specified by the direction cosines with respect to the 
x (black), y (blue), and z (green) directions. For the MFID 
simulations, this (x, y, z) coordinate system was rotated in an 
arbitrary way relative to the original simulation coordinates 
to see if the latter could be recovered in blind tests. 

[ 14 ] As suggested by Denton et al. [2010], there are two 
guidelines for picking a section of time within which to cal- 
culate the gradient directions. We look for a period of time 
during which there is a large separation in the gradient 
eigenvalues (Figure 4a) and during which there is also sta- 
bility to the gradient directions (Figures 4b-4d). These cri- 
teria are fairly well satisfied in the right half of the plot from 
t = 40 to 140 s. Using this time interval, we calculate the 
average perturbed gradient matrix (6M V ' S ) and then get the 
eigenvalue directions using this matrix. The direction 
cosines relative to the simulation (x, y, z) coordinate system 
calculated using this average matrix are the filled circles 
(dots) plotted in Figures 4b^4d, where the colors have the 
same meaning as for plotting the instantaneous direction 
cosines. Clearly the direction cosines based on the average 
matrix are consistent with the values calculated using the 
instantaneous gradient <5M Vfl during this time interval ( t = 
40 to 140 s). These directions are 2.7, 4.4, and 3.5° off from 
the presumed exact directions. (Since the actual simulation 



t(s) 


Figure 3. Maximum gradient (black), intermediate gradi- 
ent (blue), and minimum gradient (green) eigenvalues, A G _„ 
versus time in seconds for Run 1 using the matrix M vs 
based on the gradient of the unperturbed magnetic field. 
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tion Co-max, (c) intermediate gradient direction e G _ int , an d 
(d) minimum gradient direction e G _ m j n , all versus time for 
Run 1 . The gradient directions are specified by the direction 
cosines with respect to the x (black), y (blue), and z (green) 
directions. The filled circles (dots) indicate the directions 
cosines found using the perturbed gradient matrix (6M VS ) 
averaged from t = 40 to 140 s. 

was only quasi-two-dimensional, the exact simulation max- 
imum gradient, intermediate gradient, and minimum gradi- 
ent directions are not known with precision.) These errors 
are listed in Table 1. 


[ 15 ] Now we return to the question of how to determine 
the time interval for averaging (M VB ) 0 . Denton et al. [2010] 
suggested that the time interval should encompass the region 
of maximum gradient. This is the region where the recon- 
nection magnetic field changes rapidly within the “magnetic 
island”, the region of reconnected magnetic flux shown, for 
instance, in Figure 1. They showed [Denton et al, 2010, 
Table 1] that the results for the directions were relatively 
insensitive to the averaging interval. The results became 
somewhat worse if a very large averaging interval was used; 
the accuracy for small averaging intervals wasn’t explored. 
Now we attempt to find a systematic method to determine a 
good averaging interval. 

[ 16 ] We first identify the largest range of time for which 
the time-dependent maximum gradient is within half of its 
maximum value maximized with respect to time, i.e., the 
time range for which the black curve in Figure 3 is within a 
factor of 4 of its largest value (note that the eigenvalue is the 
square of the gradient). Then we examine the error in the 
minimum gradient direction (which direction is often diffi- 
cult to find; the maximum gradient direction can usually be 
found from the unmodified Shi et al. [2005] method) for 
different time intervals centered on the midpoint of this time 
interval. The range of time used, At range , can be specified by 
the ratio Af range /At half , where A/ half is the previously defined 
time range encompassing the region for which the maximum 
gradient is half of its maximum value. We use the same time 
range t = 40 to 140 s for evaluating the average perturbed 
gradient matrix ((5M VS ), from which we evaluate the gra- 
dient directions, in each case. Figure 5a shows the angle 
^(e G -min, e Y) between the minimum gradient direction 
^G— min and the presumed exact direction e Y . While there are 
possibly fortuitous choices of small averaging intervals, the 
presumed error Z(e G _ min , e Y ) becomes consistently small 
only for At range /At half > 2. The choice of —41.5 s and 34.5 s 
for the averaging interval mentioned earlier represented 
At range / A f ha if = 4. (Figure 5b is discussed in Section 3.2, the 
bold curve of Figure 5c is discussed in Section 3.3, and the 
thin curves of Figure 5c are discussed in Section 4.) 

[ 17 ] While in our experiments we know the desired 
answer, i.e., the e Y direction, this will not be known when we 
are dealing with real spacecraft data. But Figure 5a suggests 
that we can use a time interval for the average gradient that is 
large enough such that the minimum gradient direction 
converges to some stable answer. 

[is] For comparison, if we instead use the unperturbed 
matrix M v,s based on B with errors within the time interval 
t = — 10 to 10s in Figure 3, during which time the separation 
in the three eigenvalues is largest, the eigenvector directions 
from the average matrix are 1.0, 35.6, and 35.6° off from the 
exact directions. This shows that the maximum gradient 


Table 1 . Angular Difference Between Directions e G -max, e G -i n t, and e G -min Inferred by the Modified Shi et al. [2005] Method Using the 
Perturbed Gradient, and the Corresponding Approximate (for Quasi-Two-Dimensional MHD Simulations) or Exact (for Exactly Two- 
Dimensional Kinetic Simulation) Directions, ez, ex, and e Y 


Run 

Equations 

Guide Field 

Time Range (s) 

max? &z) 

A e G-int, exf 

min? & y ) 

1 

MHD 

no 

40 to 140 

2.7° 

4.4° 

3.5° 

2 

MHD 

yes 

-700 to -500 

3.9° 

3.7° 

3.1° 

3 

Kinetic 

no 

150 to 190 

2.6° 

IT 

1.9° 


a Angle between the argument vectors. 
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At 

range half 


Figure 5. Angle Z(e G _ mm , e Y ) between the minimum gra- 
dient direction ec-mm and e Y versus the normalized time 
range At range /A?haif used for averaging (M vs ) 0 as described 
in the text for (a) Run 1, (b) Run 2, and (c) Run 3 (thick 
black curve). In Figure 5c, the thin curves show variations 
of Run 3; the thin black curve is for the spacecraft passing 
through the X point with same velocity as was used for 
Run 3, the blue curve is for spacecraft velocity in the Z direc- 
tion crossing the current sheet about 20 c/w p ; downstream of 
the X point, and the green curve is for the spacecraft relative 
velocity in the Z direction passing through the X point. 

direction can be accurately determined using the unper- 
turbed matrix, but the minimum gradient direction can only 
be well determined (at least with the Shi et al. [2005] 
method) using the perturbed matrix if the calibration errors 
are large. The effect of calibration errors on the unmodified 
method will be examined more thoroughly in the Discussion 
section. As mentioned by Denton et al. [2010], comparison 
of the directions using the unperturbed and perturbed matrix 
can help to establish whether or not the calibration errors are 
significant. 

[ 19 ] Now we calculate the structure velocity V str using the 
data from Run 1, again with 0.01 nT random noise errors 
and 0.1 nT systematic calibration errors added to each 
component of the magnetic field. In order to calculate Vstr, 
we use the perturbed gradient and perturbed time derivative 
as described by (3) with the same time interval for averaging 


as was used for determining the directions. Before calculat- 
ing the gradients and time derivative of the magnetic field, 
the magnetic field values were averaged at each point in time 
over a time interval of ±5 s (11 data points) as before. 

[ 20 ] Figure 6 shows the results for the inferred velocity 
components versus time (curves) and the exact values of the 
components from the simulation (dotted horizontal lines). 
The components of the spacecraft measurements in the x, y, 
and z coordinate directions are indicated by black, blue, and 
green color, respectively. As was done by Denton et al. 
[2010], we determine the total structure velocity Vstr, not 
just the “in-plane” components. This is possible because the 
simulation was weakly three-dimensional and the gradient in 
the minimum gradient direction was large enough so that the 
component of the velocity in that direction could be well 
determined. Evidence that the total velocity is reliable is that 
the components of V str in Figure 6 are relatively stable 
(except in the vicinity of the spike in K str components at t ~ 0 
where the minimum gradient is particularly small; see 
Figure 3a). While we were able to determine all three com- 
ponents of the structure velocity in this case, it should be kept 
in mind that the component of the velocity in the minimum 
gradient direction could be inaccurate [Shi et al, 2006]. 

[ 21 ] Median values of the velocity components were 
determined between t = 35 s and t = 125 s, during which time 
the velocity components were particularly stable. The median 
values are plotted in Figure 6 as the filled circles. The mag- 
nitude of the inferred velocity is accurate to 1.6%, and the 
direction is accurate to within 5.2°. These values are listed in 
Table 2. (Comparable values could be found by using a time 
range which included the sharp peaks in Figure 6 as was done 
by Denton et al. [20 10]; using the median values excludes the 
extreme values within the peaks.) 



Figure 6. Run 1 spacecraft measurements x (black curve), 
_v (blue curve), and z (green curve) components (arbitrary 
coordinate system) of the structure velocity K str versus time, 
found from the Shi et al. [2006] method modified by using 
the perturbed gradient and time derivative of B. The filled 
circles show the median values of the curves using the data 
between t = 35 s and 125 s, while the dotted horizontal lines 
show the exact values from the simulation. 
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Table 2. Difference A(F str , V^f) Between the Magnitude of the 
Actual Structure Velocity V str and That Inferred Using VB, 
and the Angular Difference Z(e Vstn eyftr) Between e V str and evitr 
for the Cases Indicated 


Run 

Equations 

Guide 

Field 

Time 
Range (s) 

A (V sll , V™) 

^(^Vstn ^Usfr) 

1 

MHD 

no 

35 to 125 

1.6% 

5.2° 

2 

MHD 

yes 

-700 to -200 

1.5% 

0.9° 

2 

Kinetic 

no 

150 to 180 

0.08% 

1.8° 


3.2. MHD Simulation With a Guide Field 

[ 22 ] Run 2 is also quasi-two-dimensional. It differs from 
Run 1 in that it has a guide field (a component of B in the 
“out of plane” direction) of about half the magnitude of the 
reconnecting field. Again we add 0.01 nT random noise 
error and 0.1 nT systematic calibration errors to each com- 
ponent of the magnetic field. Next we subtract the average 
gradient within t = — 68 and t = 60 s to determine the per- 
turbed gradient. In order to make this choice, we examined 
the presumed error in the minimum gradient direction in a 
similar manner as was done for Run 1. The results 
corresponding to A/ range /At half = 8 are shown in Figure 5b. 
(As was mentioned in reference to Run 1, when we are 
analyzing spacecraft data we will not know the Y direction, 
so in practice, we will need to look for convergence of the 
minimum gradient direction.) Figure 7a shows the eigenva- 
lues of the maximum, intermediate, and minimum gradient 
directions for Run 2, while Figures 7b-7d show the 
corresponding direction cosines (as for Run 1). Again, we 
look for a time interval during which the eigenvalues for the 
gradient directions (Figure 7a) are well separated and the 
eigenmode directions are relatively steady. Picking the time 
period between t = —700 and —500 s, we find the compo- 
nents (black, blue, and green for x, y, and z, respectively) of 
each eigenmode direction; these are plotted as the filled 
circles in Figures 7b-7d. These directions are 3.9° 3.7° and 
3.1° off from the presumed extremal gradient directions in 
the simulation; these values are listed in Table 1. As men- 
tioned in reference to Run 1, the true directions are not 
exactly known because the simulation is only approximately 
two-dimensional. If we instead choose the time period t = 40 
to 90 s, for which the eigenvalues are better separated 
(Figure 4a), but the directions are less steady (Figures 7b-7d), 
we find that the directions are 5.8°, 7.0°, and 5.4° off from the 
presumed true directions. 

[ 23 ] If we do not add any random or calibration errors to 
the simulation data, and use the unmodified Shi et al. [2005] 
method, we get the results shown in Figure 8. In this case, it 
is not possible to use the same time interval to find the 
maximum and minimum gradient directions. Using the 
central time interval, t = — 15 to 5 s, where the maximum 
gradient eigenvalue (black curve in Figure 8a) is extremely 
large, we can find the maximum gradient direction ^G-max 
with an error of 3.4° from the presumed direction. Then 
using the interval t = 120 to 260 s, during which the mini- 
mum gradient eigenvalue Ao-min (green curve in Figure 8a) 
is much lower than the other eigenvalues, and during which 
the minimum gradient eigenmode direction <? (l -min is rela- 
tively steady (Figure 8d), we can find the minimum gradient 
direction with an approximate error of 5.7° from the pre- 
sumed direction. If we then get the intermediate gradient 


direction as the renormalized cross product e G - m in x 
^G-max, we find the presumed errors in the three directions 
to be 3.4°, 1.0°, and 4.1°, respectively, for the maximum, 
intermediate, and minimum gradient directions. 

[ 24 ] Comparison of Figures 7 and 8 shows an interesting 
fact. Subtracting off the maximum gradient near the peak of 
the maximum gradient (t ~ 0 in Figure 8a) has the effect of 
introducing a larger maximum gradient (in that direction) at 
times away from the peak (times different than / ~ 0 in 
Figure 4a; note that the black curve is almost uniformly at a 
large value). So, in addition to subtracting off constant gra- 
dients that would result from calibration errors, the modified 
Shi et al. [2005] method using the perturbed gradient trans- 
fers information about the maximum gradient to times away 
from its peak value. 

[ 25 ] Now we calculate the structure velocity V stT using the 
data from Run 2 with 0.01 nT random noise error and 0.1 nT 
systematic calibration error added to each component of the 
magnetic field. We then smooth the data at each time using a 
time interval ±50 s (101 data points) and use the perturbed 



t(s) 


Figure 7. Same as Figure 4, except Run 2 (MHD with a 
guide field). The filled circles indicate the directions cosines 
found using the perturbed gradient matrix (hM VB ) averaged 
from t = —700 s to —500 s. 
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Figure 8. Same as Figure 7, except for Run 2 with no mag- 
netic field errors added, and analyzed using the unmodified 
Shi et al. [2005] method. 

gradient and perturbed time derivative as described by (3) to 
calculate The results are shown in Figure 9 for the entire 
velocity V str , again not subtracting off the component of ^sh- 
in the minimum gradient direction. The clue that we are able 
to obtain a reliable result is again that all 3 velocity com- 
ponents in the plots are steady (Figure 9). Figure 9 shows the 
results for the inferred velocity components versus time 
(curves) and the exact values of the components from the 
simulation (dotted horizontal lines). The median values of 
the velocity components were determined between t= — 700 
and —200 s, and are plotted in Figure 9 as the filled circles. 
The magnitude of the inferred velocity is accurate to 1.5%, 
and the direction is accurate to within 0.9°. These values are 
listed in Table 2. 

[ 26 ] Figure 10 shows the inferred structure velocity Vstr 
using the same method, but with a smoothing interval of 
±5 s instead. The results are similar, but more noisy. 

3.3. Full Dynamics Kinetic Simulation 

[ 27 ] Now we consider results from an array of virtual 
spacecraft flying through a two-dimensional full dynamics 
kinetic simulation of the magnetotail. This simulation was 


described by Shay et al. [2007] and Drake et al. [2009]. In 
this case, referred to as Run 3, the simulation is truly two- 
dimensional. While there is no constant guide field, there is a 
significant spatially varying Flail magnetic field in the out of 
plane direction, with magnitude about half that of the 
reconnection magnetic field immediately upstream of the 
current sheet. The kinetic simulation used normalized units; 
for the purposes of comparison and more easily relating to the 
real situation in the magnetotail, we convert the normalized 
units to real units by multiplying values of B by 15 nT, values 
of time by (1/1.5) s, values of velocity by 1500 km/s, and 
values of distance by 1000 km. If the fields were exactly 
known, it would be possible to find the direction along which 
there was exactly zero gradient. The simulation, however, 
has numerical noise due to the particle distribution. We 
estimate the amplitude of this numerical noise level to be at 
least 0.005 nT. In order to bring the noise level up to 0.01 nT, 
we add additional time-dependent random noise (time- 

dependent) of 0.009 nT (from \J (0.01 ) 2 — (0.005) 2 ) to each 
component of the magnetic field. In addition, we also add 
random time independent calibration errors to the compo- 
nents of the magnetic field at the level of 0.1 nT. The data 
were sampled at a time resolution of 0.1 s, and the effect of 
noise errors could be eliminated with a smaller averaging 
interval; we used ±3 s. 

[ 28 ] The interval used to determine the average gradient, 
t = 70 to 120 s, was chosen in the same way as for Runs 
1 and 2. The error in the minimum gradient direction is 
shown versus the averaging interval as the thick black 
curve in Figure 5c. We chose At range /At half = 4. 

[ 29 ] Figure 11a shows the eigenvalues of the maximum, 
intermediate, and minimum gradient directions for Run 3, 
while Figures 1 la-1 Id show the corresponding direction 
cosines using the modified Shi et al. [2005] method. 

[ 30 ] The eigenvalues are well separated and the directions 
are stable in the early and late stages of the simulation. Using 
the time period between 150 and 190 s, we find directions 
that are 2.6°, 2.7°, and 1.9° off from the exact directions 



Figure 9. Like Figure 6, except for Run 2. The filled cir- 
cles show the median values of the curves using the data 
from t = —700 s to —200 s. 
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Figure 10. Same as Figure 9, but with the smoothing inter- 
val for the magnetic field reduced to ±5 s. 

(because this simulation is precisely two-dimensional); these 
values are listed in Table 1. In this case, the e x , e y , and e z 
spacecraft measurements directions are the same as the 
simulation directions e x , e Y , and e z , respectively (the coor- 
dinate axes have not been rotated). Clearly, these directions 
are very well determined by the modified Shi et al. [2005] 
method. 

[ 31 ] If we calculate the eigenvalues and directions using 
the unmodified Shi et al. [2005] method instead, but still 
including random and calibration errors, the directions (not 
shown) are not very steady for any time period. The most 
stable directions occur near the peak in the gradient, but the 
intermediate and minimum gradient directions found using 
that time interval are off from the actual directions by about 
88°. The maximum gradient direction is accurate to within 
1.5°. This again shows that calibration errors at a level of 
0.1 nT severely corrupt information about the intermediate 
and minimum gradient directions (which have relatively 
small contributions to the total gradient, compared to the 
maximum gradient direction), but do not significantly con- 
taminate information about the maximum gradient direction. 

[ 32 ] Figure 12 shows the same data as in Figure 11, 
except that here the unmodified Shi et al. [2005] method is 
used for the Run 3 data without errors. The filled circles in 
Figures 12b— 12d show the directions plotted in Figure 11 
using the modified method; these directions are within 
about 3° of the exact directions (Table 1), as was discussed 
previously. This figure also shows some of the features we 
have discussed previously. The maximum eigenvalue direc- 
tion is steady and pointing in the e z ( =e z ) direction (the green 
curve in Figure 12b is close to unity) within the central time 
region from about f = 70 to 125 s, where the maximum gra- 
dient is large ( the black curve in Figure 12a is well above the 
blue and green curves). The correct direction of the minimum 
gradient is the e v (=e y ) direction, as indicated by the blue 
curve in Figure 12d being near plus or minus 1 (there is an 
arbitrary 180° ambiguity in the direction). The minimum 
gradient direction in Figure 12d is seen to be accurate in the 
early (before 50 s) and late (after 125 s) times. These intervals 
correspond to the times during which the minimum gradient 


eigenvalue (green curve in Figure 12a) is significantly lower 
than the other eigenvalues (blue and black curves in 
Figure 12a). During the period of large maximum gradient, 
the minimum gradient direction from the unmodified Shi 
et al. [2005] method is not steady, indicating that the best 
time periods for determining the minimum gradient direction 
are at early or late times. 

[ 33 ] Now we calculate the structure velocity for Run 3 
using the modified Shi et al. [2006] method, again including 
calibration and noise errors. Since Run 3 is truly two- 
dimensional, the gradient in the out of plane direction is 
exactly zero. Therefore the out of plane component of F str is 
arbitrary, and cannot be found from (3), except as deter- 
mined by the noise. Figure 13 shows the results for ^str 
using (3) to determine only the components in the maximum 
and intermediate gradient directions. As described earlier, 
we average the gradients and dB/dt in the interval 1 = 70 to 
120 s and also smooth the fields over an interval of 
approximately ±30 s (increasing this averaging interval does 
not significantly improve the results). Figure 13 shows that 



Figure 11. Same as Figure 4, except for Run 3 (the full 
dynamics kinetic simulation). The filled circles indicate the 
directions cosines found using the perturbed gradient matrix 
(bM VB ) averaged from t = 150 s to 190 s. 
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Figure 12. Same as Figure 11, except that the unmodified 
Shi et al. [2005] method (using the total gradient of B) is used 
for the Run 3 data with no magnetic field errors added. The 
filled circles in Figures 12b— 12d show the directions from 
Figure 11, i.e., from the modified method, again for t = 
150 s to 190 s. 


the in plane components of Vstr (black and green curves) are 
well determined. The median values in the interval 150 to 
180 s are plotted as filled circles, and agree well with the 
exact values (horizontal lines), with an error of 0.08% in 
magnitude and 1.8° in direction. These values are listed in 
Table 2. Note that V strj , in Figure 13 (blue curve) is not 
exactly zero because the minimum gradient direction used to 
subtract off the minimum variance component of Vstr is the 
one that we calculated and not the exact Y direction (here the 
y direction). 

[ 34 ] Figure 14 shows the results if we try to determine all 
three components of Var. The x and z (in plane) components 
of F str determined from (3) (black and green curves, with 
median values shown for the time interval t = 150 to 180 s as 
filled circles) again agree well with the exact values (dotted 
horizontal lines) but the out of plane component K str v (blue 
curve in Figure 14) is not accurate, despite the fact that it 
appears to be steady from t = 150 s to 180 s (though not at 


earlier times). It should have the same value as the z com- 
ponent (green curve). This indicates that caution should be 
used when interpreting the minimum gradient component of 
the structure velocity, at least using the modified Shi et al. 
[2006] method. (For the unmodified method, the y compo- 
nent of Vstr varies wildly at all times (not shown), clearly 
indicating that it is not meaningful. ) 

4. Discussion and Summary 

[ 35 ] We have shown that the Shi et al. [2005, 2006] 
methods can be successfully used to determine the recon- 
nection geometry and structure velocity Vstr relative to the 
spacecraft, even when calibration and noise errors are pres- 
ent, provided that the perturbed gradient is used as described 
in Section 1. This was shown in Section 3.1 for an MHD 
simulation without a guide field similar to that examined by 
Denton et al. [2010], in Section 3.2 for an MHD simulation 
with a guide field, and in Section 3.3 for a full dynamics 
kinetic simulation that did not have a constant guide field, 
but did have a Hall magnetic field in the out of plane 
direction. The agreement between the inferred directions and 
the approximate true directions (for quasi-two-dimensional 
simulations) or the exact directions (for the two-dimensional 
simulation) was good (all within 5°, as shown in Table 1). 
The structure velocity could also be well determined (mag- 
nitude within 2% and direction within 5°, as shown in 
Table 2). 

[36] The maximum gradient and the associated large 
eigenvalue of the gradient matrix occurs generally at the 
center of the current sheet crossing, whereas the minimum 
eigenvalue is most distinct from the intermediate eigen- 
value away from the central current sheet crossing (Figures 3 
and 12). Thus if the unmodified Shi et al. [2005] method is 
used to determine the reconnection geometry, our results 
suggest that the maximum gradient direction may be best 
determined near the central current sheet, while the minimum 



t(s) 


Figure 13. Like Figure 6, except for Runs 3. The compo- 
nents of the velocity are determined only in the maximum- 
intermediate gradient plane. The filled circles show the 
median values between t = 150 s and 180 s. 
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Figure 14. Like Figure 13 for Run 3, except that now all 
the components of the velocity are determined. The filled 
circles show the median values between t— 150 s and 180 s. 

gradient direction may be best determined away from the 
central current sheet. 

[ 37 ] To use the modified Shi et al. [2005, 2006] methods, 
we need to subtract the average gradient from a particular 
time interval. In agreement with results of Denton et al. 
[2010], we find that a good choice for the subtracted term 
(M vs ) 0 in (2) is the average gradient near the central current 
sheet crossing. Examination of Figures 3 and 12 shows that 
the maximum eigenvalue tends to be best separated in 
magnitude from the other eigenvalues (intermediate and 
minimum) in this region, and that the minimum and inter- 
mediate eigenvalues are often least well separated. This 
means that, at the central current sheet crossing, the recon- 
nection structure tends to be more one dimensional than at 
other locations. After subtracting the gradient in the central 
current sheet, information about that gradient direction is 
transferred to regions far from the central current sheet 
(Figures 4a, 7a, and 1 la), so that an interval away from the 
central current sheet can then be used to determine all three 
directions (Section 3.1). 

[38] The structure velocity F str can also be determined 
using the modified Shi et al. [2006] method. If the structure 
is three-dimensional, we are able to determine all three 
components of the velocity (Figures 6 and 9), but if the 
minimum gradient is very small (such as occurs exactly for a 
true two-dimensional system), the minimum gradient com- 
ponent of the structure velocity will not be well determined. 
In that case we can only get the “in-plane” components of 
V s « (Figures 13 and 14). Noise errors lead to noise in the 
inferred values of V stT (Figure 10 versus Figure 9). The time 
resolution of the magnetic field instrument is greater for 
CLUSTER (1/67 s 1 time resolution [Balogh et al., 1997]) 
and even greater for MMS (0.01 s resolution for the fluxgate 
instrument), than what we have assumed here (between 0.1 
and 1 s). We therefore conclude that digital noise errors are 
unlikely to be a problem. 

[ 39 ] In all the cases we have studied up to this point, the 
virtual spacecraft passed to the side of the X point as sche- 
matically illustrated in Figure 1. For the full dynamics 


kinetic simulation, we tried some variations of the spacecraft 
trajectory. For Run 3, the spacecraft trajectory was at an 
angle of about 27° to the current sheet and passed about 
20 du p i (assumed to be ~20,000 km) in the outflow 
downstream of the X point. The thin black curve in 
Figure 5c shows the error as a function of the averaging 
interval if the spacecraft trajectory is at the same angle to the 
current sheet but now shifted in the X direction so as to pass 
through the X point. The blue curve in Figure 5c shows the 
error as a function of the averaging interval if the spacecraft 
trajectory is normal to the current sheet (in the Z direction) 
and crosses the current sheet about 20 du pi downstream of 
the X point. The green curve in Figure 5c shows the error as 
a function of the averaging interval if the spacecraft trajec- 
tory is normal to the current sheet (in the Z direction) and 
crosses the current sheet through the X point. Though some 
of the errors in Figure 5c (^10°) are somewhat larger than 
those in Table 2 (<5°), it appears that the modified Shi et al. 
[2005] method can be used for all these cases as long as the 
time interval for averaging the gradient matrix is sufficiently 
large around the region of peak gradient. We need to note 
one complication. Note that the blue curve in Figure 5c does 
not extend past values of Af range /At half = 4. For greater 
values, the minimum gradient direction is not clearly steady 
in any time interval, and in that case the results will depend 
heavily on the time interval used. 

[ 40 ] We found that one needs to examine carefully the 
effect of different values of Af range /A? half , especially for the 
kinetic simulation. To do this we used the fact that we knew 
the right answer. In order to test again how we would do in a 
truly blind test, we generated a totally new set of data from 
the simulation used for Run 3. In this case, the spacecraft 
trajectory was at an angle of 45° to the current sheet passing 
10 c/iOpi downstream of the X point. The coordinates were 
then rotated using random Euler angles that were unknown 
to us. We then chose a time interval during which the three 
directions were relatively constant, and found the maximum, 
intermediate, and minimum gradient directions for different 
values of At range /At half . We looked for convergence of the 
directions. The directions for At range /At half = 8 were within 
0.1° of the directions using At range /At half = 32, while the 
directions for At range /At half = 4 varied by about 5° from 
those using At range /At half = 32. Based on this, the directions 
were determined using At range /At half = 8. After rotating the 
resulting directions back to the original coordinate system, 
the directions for maximum, intermediate, and minimum 
gradient direction were found to be 9.6°, 9.7°, and 1.5° off 
from the correct directions. These results were obtained 
using the modified Shi et al. [2005] method from the per- 
turbed gradient. We also calculated the maximum gradient 
direction using the unmodified method (without subtracting 
the average gradient). In that case, we found the maximum 
gradient off by only 0.9° from the correct direction. 

[ 41 ] Figure 15 shows the median error for 10 trials of 
Run 3 with different random calibration errors of varying 
magnitude and different noise errors of the same magnitude 
(0.01 nT) for the maximum gradient direction e G - max (black 
curve) and minimum gradient direction e G -mm (blue curve) 
versus the ratio of the calibration error dB cal to the recon- 
nection magnetic field B rec expressed as a percentage. For 
each trial, a random distribution of dB was assumed for each 
component of the magnetic field. The horizontal axis shows 
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Figure 15. Median error for 10 trials with different random 
calibration of differing amplitude and noise errors of the same 
amplitude for the maximum gradient direction <Y,_ max (black 
curve) and minimum gradient direction e G -min (blue curve) 
versus the calibration error dB cal . This error is expressed as 
a percentage of the reconnection magnetic field B rcc . 

the root mean squared value dB ca i for each component of the 
assumed distribution. Based on Figure 15, the unmodified 
Shi et al. [2005] method can determine the maximum gradi- 
ent direction to within 3° for calibration errors as large as 3% 
of the reconnection magnetic field, while in order to find the 
minimum gradient direction to the same accuracy, the cali- 
bration errors must be no larger than 0.03% of the recon- 
nection magnetic field. For reference, the calibration errors 
assumed here are 0.7% of the reconnection magnetic field. 

[ 42 ] Based on these numbers, it’s likely that the unmodi- 
fied method can be used to determine the maximum gradi- 
ent direction, but that the modified method must be used 
to detennine the minimum gradient direction. The results 
described above for the blind test suggest that the maximum 
gradient direction might be determined more accurately by 
the unmodified method. The method is also less complex, 
not depending on a choice of time interval for averaging the 
gradient matrix. Combining the maximum gradient direction 
from the unmodified method with the minimum gradient 
direction from the modified method, all three angles would 
be determined within a couple degrees. Another possibility 
would be to combine the maximum gradient direction from 
the unmodified method with the maximum variance direction 
for the magnetic field (yielding the intennediate gradient 
direction), and from those two directions get the minimum 
gradient direction. 

[ 43 ] It should be kept in mind when using the modified Shi 
et al. [2005, 2006] method that subtraction of the gradient 
(M vb ) 0 averaged over the interval Af range (in order to cancel 
out the effect of constant (in time) calibration errors) results 
in loss of real information about the average gradient of the 
real field. In a situation with constant (in time) gradients, the 
modification would not be likely to give satisfactory results. 
Comparison between results from the unmodified Shi et al. 
[2005, 2006] methods and the modified methods using the 


perturbed gradient and time derivative of B may give infor- 
mation on the degree to which calibration errors contaminate 
the gradients. Ideally, a way will be found to reduce the 
calibration errors sufficiently as mentioned in Section 2 so 
that the unmodified and modified methods using the per- 
turbed fields yield the same results. But our estimates indi- 
cate that if these calibrations are not reduced, results from 
the unmodified method may not be satisfactory, especially 
for the minimum gradient direction (Figure 15). 

[ 44 ] The effects of explicit time dependence still need to be 
investigated. It’s likely that explicit time dependence will 
make the procedures we have demonstrated here more diffi- 
cult. Directions and velocities may be changing with respect 
to time, in which case they will have to be determined within 
smaller time intervals. In that case, we would not be able to 
use extremely large values of At range /At half . Presumably, as 
At ra nge/ A t ha if is increased, the directions would converge 
toward some values, and then diverge as the timescale over 
which the data is collected becomes large. Also, large 
amplitude time-dependent effects in a 3-D kinetic system, 
associated with physical effects such as waves or non- 
uniformities in the inflow to the reconnection site, may be 
larger than what we have in our simulations, and these could 
have an adverse effect on our ability to obtain accurate 
information from the Shi et al. [2005, 2006] methods and 
their extension to perturbed fields [Denton et al., 2010], 
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